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A new method is introduced to create artificial time se- 
quences that fulfil given constraints but are random other- 
wise. Constraints are usually derived from a measured signal 
for which surrogate data are to be generated. They are ful- 
filled by minimizing a suitable cost function using simulated 
annealing. A wide variety of structures can be imposed on the 
surrogate series, including multivariate, nonlinear, and non- 
stationary properties. When the linear correlation structure 
is to be preserved, the new approach avoids certain artifacts 
generated by Fourier-based randomization schemes. PACS: 
05.45.-Hb 

Randomization of data and Monte Carlo resampling of 
probability distributions is a common technique in statis- 
tics [Q. In the context of nonlinear time series analysis 
it has been discussed by several authors and is usually 
referred to as the method of surrogate data Q|. A null 
hypothesis for the nature of a time series can be tested 
by comparing the value of an observable 7 obtained us- 
ing the data with values obtained using a collection of 
surrogate time series representing the null hypothesis. 
All but the simplest null assumptions allow for certain 
structures, for example linear serial correlations. There 
are two distinct ways to implement such structures when 
creating surrogate series. Traditional bootstrap methods 
use explicit model equations that have to be extracted 
from the data. This typical realizations approach can 
be very powerful for the computation of confidence in- 
tervals, provided the model equations can be extracted 
successfully. As discussed in Ref. the alternative ap- 
proach of constrained realizations is more suitable for the 
purpose of hypothesis testing. It avoids the fitting of 
model equations by directly imposing the desired struc- 
tures onto the randomized time series. However, the 
choice of possible null hypothesis has so far been lim- 
ited by the difficulty of imposing arbitrary structures on 
otherwise random sequences. Algorithms exist mainly 
for the following cases. (1) The null hypothesis of in- 
dependent random numbers from a fixed but unknown 
distribution can be tested against permutations without 
repetition of the data since these conserve the sample dis- 
tribution exactly. (2) The case of Gaussian noise with ar- 
bitrary linear correlations leads to the Fourier transform 
method. The Fourier transform of the data is multiplied 
by random phases and then transformed back, conserving 
the sample periodogram. (See Ref. [Q for the multivari- 
ate case.) (3) Surrogates with a given distribution and 
given linear correlations are needed for the null hypothe- 
sis of a monotonically rescaled Gaussian linear stochastic 
process. This is approximately achieved by the amplitude 



adjusted Fourier transform (AAFT) algorithm |Q| and the 
more accurate iterative method proposed in Ref. Q. 

This paper will introduce a general method for gen- 
erating random time sequences subject to quite general 
constraints. Any null hypothesis that leads to a complete 
set of observables can thus be tested for. All the above 
cases can be dealt with (often with higher accuracy) , but 
also multivariate, nonstationary, nonlinear or other con- 
straints can be implemented. In all the applications in 
this paper, the single time probability distribution will 
be one of the constraints, leading to the requirement that 
the randomized sequence is a permutation of a fixed col- 
lection of values. All other constraints, for example part 
or all of the lags of the autocorrelation function, will be 
formulated in terms of a cost function which is then min- 
imized among all possible permutations by the method 
of simulated annealing. 

After giving the actual randomization scheme I will 
discuss some major applications. We will show that the 
algorithm yields a more accurate nonlinearity test and 
avoids known artifacts that are introduced by end effects 
with ordinary, Fourier-based surrogates . We will also 
give examples with more general null hypothesis than 
that of a rescaled stationary linear stochastic process. 
For these examples, previous methods could not provide 
appropriate surrogates. 

The algorithm is conceptually very simple: 

1. Specify constraints Ci({in}) = in terms of a 
cost function i5({i„}), constructed to have a global 
minimum when the constraint is fulfilled. 

2. Minimize E{{xn}) among all permutations of 
a time series {xn} by simulated annealing. Config- 
urations are updated by exchanging pairs in {x„}. 

Examples of its use will be given below. 

The simulated annealing method is particularly use- 
ful for combinatorial minimization with false minima. It 
goes back to Metropolis et al. [|j , and is thoroughly dis- 
cussed in the literature . Essentially, the cost function 
is interpreted as an energy in a thermodynamic system. 
At some finite "temperature" T, system configurations 
are visited consecutively with a probability according to 
the Boltzmann distribution e^^^'^ of the canonical en- 
semble. This is achieved by accepting changes of the 
configuration with a probability p — \ li the energy is 
decreased (A£' < 0) and p ~ ^-^^1^ if the energy is in- 
creased, (AE' > 0). The temperature is decreased slowly, 
thereby "annealing" the system to the ground state of 
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minimal "energy" , that is, the minimum of the cost func- 
tion. In the hmit T — > 0, all ground state configura- 
tions can be reached with equal probability. Although 
some general rigorous convergence results are available, 
in practical applications of simulated annealing some 
problem-specific choices have to be made. In particular, 
apart from the cost function itself, one has to specify a 
method of updating the configurations and a schedule for 
lowering the temperature. A way to efficiently reach all 
permutations by small individual changes is by exchang- 
ing randomly chosen (not necessarily close-by) pairs. In 
many cases, an exchange of two points is reflected in a 
rather simple update of the cost function. This is impor- 
tant for speed of computation. Many cooling schemes 
have been discussed in the literature In this work, 
the temperature is multiplied by a at each cooling step. 
Cooling is done if either the number of successful up- 
dates since the last cooling exceeds Ngucc, or the total 
number of configurations visited during this cooling step 
exceeds iVtotai- It is difficult to give general rules on how 
to choose a, Ngucc, and A'^totai- Slow cooling is necessary 
if the desired accuracy of the constraint is high. It seems 
reasonable to increase Ngucc and A'totai with the system 
size, but also with the number of constraints incorpo- 
rated in the cost function. Generally, one can choose a 
tolerance for the constraints, start with rather fast cool- 
ing and repeat the analysis with a slower cooling rate 
if the accuracy has not been met. Other more sophisti- 
cated cooling schemes may be suitable depending on the 
specific situation. The reader is referred to the standard 
literature 

Let us first demonstrate that the algorithm yields more 
accurate results than previous methods for the most 
prominent application of surrogate data, which is sta- 
tistical testing for nonlinearity in a time series. Con- 
sider the null hypothesis that there is a sequence {j/n} 
that has been generated by a Gaussian linear stochastic 
process. As the only allowed kind of nonlinearity, the 
actual data {a;„} consists of observations of {j/n} made 
through a monotone instantaneous measurement func- 
tion: Xn — fiUn)- As discussed e.g. in Ref. [||, the 
corresponding Monte Carlo sample has to be constrained 
to have (i) the same single time probability distribution 
and (ii) the same sample auto-covariance function jl^ 
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for all lags r = 0, . . . , TV — 1. (Zero mean has been im- 
posed for simplicity of notation.) 

In the actual test, a nonlinear observable 7 is com- 
puted for the data and a collection of surrogate data 
sets. (See Ref. 11 1 for a comparison of the performance 
of different statistics 7.) The null hypothesis will be re- 
jected if the result 70 obtained for the data is incompati- 
ble with the probability distribution of 7 estimated from 
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scramble 
AAFT 
iterative 
annealing 



0.8 
0.9 
0.98 



0.01s 
2s 
2m 

25m 
lOh 



0.82 ± 0.02 
0.08 ± 0.02 
0.03 ± 0.01 
0.0055 
0.0009 
0.0003 



TABLE I. Residual deviation from the desired auto-co- 
variance function for different methods of randomizing time 
series. 



the surrogates. Note that although even different real- 
izations of the same process will have the same sample 
auto-covariance function only up to statistical fluctua- 
tions, it is essential that the surrogates are constrained 
to C (t)'^'^'^*'^'^ as accurately as possible-since almost ev- 
ery discriminating statistic 7 will depend on C(t), we are 
otherwise likely to introduce a bias and possibly spurious 
rejections of the null hypothesis. See also the discussion 
in Ref. ||. 

Previous attempts to implement the above constraints 
have only been partially successful. In the scheme intro- 
duced here, property (i) is easily implemented by con- 
sidering as candidates for randomized series all permu- 
tations of the measured time sequence {yn\- Require- 
ment (ii) can be achieved by finding a permutation of 
{Vn} which, within the desired accuracy, minimizes a cost 
function like the following fl^j: 
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Provided the annealing scheme is brought to convergence 
with high accuracy, the known artifacts that remain with 
previous approaches can be avoided. 

As discussed in js), the original (AAFT) algorithm ||] 
can show a bias towards a flat spectrum for short se- 
quences. The iterative scheme proposed in [|] removes 
this bias to a satisfactory approximation for practical 
work. Let us however compare the accuracy of the pre- 
viously proposed schemes to the present algorithm. For 
comparability, a cost function is chosen with respect to 
the time periodic sample auto-covariance function 



n=0 



(3) 



which corresponds to the Fourier spectrum through 
the Wiener-Khinchin theorem. 



max^^o \C'p{t) - Cj 
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Minimizing Ep°°^ 
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covariance Cp"^'^^'^' measured on the data 



(r)| will reproduce the auto- 
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of length N = 1000 are generated by an autoregres- 
sivc model but measured using a nonlinear measurement 
function: Xn — Vn^ Un = 0.9?/„_i -I- 77„. The residual 
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maximal deviations of the auto-covariances of the time 
series and surrogate sets were determined for (i) random 
permutations of the data, (ii) usual AAFT surrogates , 
(iii) surrogates created with the iterative scheme given in 
Ref. 1^ and (iv) outcomes of the annealing procedure for 
different cooling protocols. Note that with slower cool- 
ing, arbitrarily high accuracy can be reached in principle. 
Averages over 20 realizations were determined for cases 
(i) to (iii). The iterative scheme (iii) was repeated until 
a fixed point was reached which was the case after about 
200 iterations. Table|| summarizes the results. Computa- 
tion time time on a DEC alpha workstation at 400MHz 
clock rate are given only for relative comparison. The 
price for the superior accuracy of the annealing scheme 
is its much higher computational cost. 

As mentioned earlier, all previous randomization 
schemes [||j5) make use of the Fourier transform in order 
to achieve the desired linear correlation structure. Note, 
however, that two sequences with the same Fourier am- 
plitudes do not quite have the same auto-covariance func- 
tion C(t), eq. (|l]). The Wiener-Khinchin theorem only 
says that the periodic sample auto-covariance function 
Cp{t), eq. (H), will be the same. This amounts to assum- 
ing that the measured time series is exactly one period of 
an infinite periodic signal, which is of course not what we 
believe to be the case. The artifact generated by this flaw 
of previous algorithms has been discussed e.g. in Ref. 
The periodically extended sequence may undergo a phase 
slip or even a finite jump at n = iV. The surrogate se- 
ries will have the power contained in that slip spread out 
over the whole observation time, leading to additional 
high frequency content. Although spurious results can be 
partially suppressed by selecting a segment of the data 
that approximately returns to the initial value, it is de- 
sirable to preserve the auto-covariance function C(t) in 
eq. (0) rather than Cp(t) in eq. (||). With the annealing 
scheme proposed in this paper, this can be easily done 
by choosing an appropriate cost function. 

As an illustration, consider a particular autoregressive 
process of order two, Xn = l.3xn-i — 0.31x„_2 -|- 
Since it is almost unstable, short realizations often show 
a large difference between the first and the last point. 
Periodic continuation turns this difference into a large 
step with broad frequency content. For a realization 
of 160 points we found that for a Fourier-based surro- 
gate (method in Ref. [||, same periodic auto-covariance 
function Cp(t)), the sample autocorrelation C(1)/C(0) 
was reduced from 0.92 to 0.85. Consequently, the power 
in the first differences is increased by a factor of two 
and short term predictability is strongly reduced. This 
can lead to spurious rejections of the null hypothesis of 
a linear process. A sequence obtained by minimizing 
E^°°'> = max^Jgi |C(t) - C'^'^'''''\t)\/t yielded the cor- 
rect value of C(l) within 2 x 10^''. 

Apart from its potential for greater accuracy, the most 
striking feature of the new scheme is its generality and 
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FIG. 1. Simultaneous measurements of breath and heart 
rates [|l^ , upper and middle traces. Lower trace: a surrogate 
heart rate series preserving the autocorrelation structure and 
the cross-correlation to the fixed breath rate series, as well as 
a gap in the data. Auto- and cross-correlation together seems 
to explain some, but not all of the structure present in the 
heart rate series. 

flexibility. This point will be demonstrated in the follow- 
ing examples which are by no means exhaustive. Note 
that none of the examples below could be studied with 
previous surrogate data schemes. Let us flrst study a 
multivariate example, a simultaneous recording of the 
breath rate and the instantaneous heart rate of a hu- 
man subject during sleep. (Data set B of the Santa Fe 
Institute time series contest in 1991 samples 1800- 
4350.) Regarding the heart rate recording on its own, 
one easily detects nonlinearity, in particular an asymme- 
try under time reversal. An interesting question however 
is, how much of this structure can be explained by linear 
dependence on the breath rate, the breath rate also be- 
ing non-time-reversible. In order to answer this question, 
one has to make surrogates that have the same auto- 
correlation structure but also the same cross-correlation 
with respect to the fixed input signal, the breath rate. 
(Here the breath rate data is not randomized, which is of 
course also possible within this framework.) Accordingly, 
a constraint is formulated involving all lags up to 500 of 
the auto-covariance and the cross-covariance {Cxy) func- 
tions. The cost bmction is taken to be max^™g |C(t) — 

other choices are possible. Further suppose that during 
one minute the equipment spuriously recorded a constant 
value. In order not to interpret this artifact as struc- 
ture, the same artifact is generated in the surrogates, 
simply by excluding these data points from the permu- 
tation scheme. 

Figure |l| shows the measured breath rate (upper trace) 
and instantaneous heart rate (middle trace). The lower 
trace shows a surrogate conserving both, auto- and cross- 
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FIG. 2. A realization of a linear process with time de- 
pendent variance (upper), a usual A AFT surrogate (middle), 
and a surrogates with the same autocorrelations and the same 
running variance as the original series. 

correlations. The cooling rate was a — 0.95, iVgucc — 
10000, A^totai = 3 X lO'^. None of the auto- and cross- 
covariances differed from the goal by more than 5 x 10^^ 
in units of the variance of the data after 3h of anneal- 
ing. (DEC alpha workstation at 400MHz clock rate.) 
The visual impression from Fig. ^ is that while the linear 
cross-correlation with the breath rate explains the cyclic 
structure of the heart rate data, other features remain 
unexplained. In particular, the surrogates don't show 
the asymmetry under time seen in the data. Possible 
explanations of the remaining structure include artifacts 
due to the peculiar way of deriving heart rate from inter- 
beat intervals, nonlinear coupling to the breath activity, 
nonlinearity in the cardiac system, and others. 

Let us finally give a nonstationary example, an AR(2) 
process with periodically modulated variance: Xn = 
1.6a;„_i — 0.8a::„_2 + &n??n with 6„ = 1 + sin^ 27m/1000. 
In Fig. H a realization {N = 2000) is shown together with 
two surrogate series. The first (middle trace) has been 
generated by the AAFT algorithm, the second (lower 
trace) has been generated by the annealing scheme to 
preserve the first 100 lags of C(t) but also the running 
variance in blocks of length 200, overlapping by 100. 

In this paper it has been demonstrated that ran- 
domization under a wide variety of constraints can be 
achieved with a permutation scheme that minimizes a 
suitable cost function using simulated annealing. The 
approach is very general. Constraints are not restricted 
to linear correlations. Multivariate, nonlinear, but also 
time dependent, nonstationary properties can be easily 
implemented. A wider range of examples will be studied 
elsewhere Q. 

Resampling with constraints is the method of choice 
for hypothesis testing, where it is preferable to paramet- 
ric bootstrap methods. Although a general, nonpara- 



metric resampling scheme has been introduced in this 
paper, care has to be taken when similar ideas are to 
be exploited for the determination of error bounds. The 
variance of statistical estimators usually depends on the 
constraints imposed. To which extent reliable error dis- 
tributions can be obtained by selecting a minimal set of 
constraints and using resampling with replacement will 
be a subject of future work. 
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